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We present a new numerical method for the construction of quasiequilibrium models of black hole- 
neutron star binaries. We solve the constraint equations of general relativity, decomposed in the 
conformal thin-sandwich formalism, together with the Euler equation for the neutron star matter. 
We take the system to be stationary in a corotating frame and thereby assume the presence of a 
helical Killing vector. We solve these coupled equations in the background metric of a Kerr-Schild 
black hole, which accounts for the neutron star's black hole companion. In this paper we adopt a 
polytropic equation of state for the neutron star matter and assume large black hole-to-neutron star 
mass ratios. These simplifications allow us to focus on the construction of quasiequilibrium neutron 
star models in the presence of strong-field, black hole companions. We summarize the results of 
several code tests, compare with Newtonian models, and locate the onset of tidal disruption in a 
fully relativistic framework. 
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I. INTRODUCTION 



Black hole-neutron star (hereafter BHNS) binaries 
have attracted considerable astrophysical interest re- 
cently. For example, they are candidates for the central 
engines of short-period, gamma-ray bursts pj . They are 
also very promising sources of gravitational radiation for 
the new generation of ground-based gravitational wave 
detectors, including the Laser Interferometer Gravita- 
tional Wave Observatory (LIGO), and the future space- 
based gravitational wave detector, the Laser Interferom- 
eter Space Antenna (LISA) . Observation of a BHNS bi- 
nary - in particular, the tidal disruption of a neutron 
star by a black hole companion - will provide spectacu- 
lar information on the gravitational fields in relativistic 
objects and the nuclear physics governing neutron star 
matter (e.g. 0). However, BHNS binaries may be among 
the most challenging objects in relativistic astrophysics 
to model numerically. Any realistic treatment must deal 
simultaneously with the difficulties associated with the 
spacetime singularity inside the black hole and the com- 
plexities associated with relativistic hydrodynamics. 

The inspiral of BHNS binaries is driven by the emis- 
sion of gravitational radiation, which extracts both en- 
ergy and angular momentum from the binary and also 
tends to circularize the orbit. As long as the binary sep- 
aration is large, the inspiral is very slow and radial ve- 
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locities are very small compared to the orbital velocities 
once the orbit has been circularized. During this phase, 
the inspiral can be modeled as a sequence of quasicircu- 
lar binary orbits. The quasiadiabatic inspiral phase ends 
when dynamical effects become important. For BHNS 
binaries, there are at least two, very different, possible 
termination points: the binary may reach an innermost 
stable circular orbit (ISCO), at which point the orbit 
becomes unstable and the black hole and neutron star 
plunge rapidly towards each other, or, prior to reaching 
the ISCO, the neutron star might be tidally distrupted 
by the gravitational field of the black hole. Predicting 
which one of these scenarios actually happens, and the 
binary separation at which it occurs, depends on a num- 
ber of factors characterizing the system, including the 
binary mass ratio. 

The importance of tidal effects can be estimated, in a 
Newtonian framework, by comparing the tidal force on a 
test mass m on the surface of the neutron star 

^jtiMbh-Rns 
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with the gravitational force exerted by the star 
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Here Mbh and Mns are the masses of the black hole and 
the neutron star, i?NS, is the neutron star radius and s 
is the binary separation. Equating the two forces yields 
an approximate critical "tidal" separation 

i? NS \m ns J ' [ ' 

or, in gravitational units with G = c = 1, 
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within which tidal forces begin to dominate and may lead 
to tidal disruption. 

Neutron stars typically have compaction ratios of 
-Rns/^ns ~ 5. That means that for neutron stars or- 
biting supermassive black holes, -Mbh 3> -WnSi so that 
the tidal separation is well inside the ISCO, which re- 
sides at about s ~ 6Mbh (the value for a point-mass 
orbiting a Schwarzschild black hole). It is therefore ap- 
propriate to neglect any effects of tidal distortion or in- 
ternal structure in calculations of neutron stars inspiral 
onto supermassive black holes. Such binaries may be a 
promising sources of low-frequency gravitational radia- 
tion for detection by LISA. 

For neutron stars orbiting stellar mass black holes, 
however, the internal structure of the neutron star may 
become important before it reaches the ISCO. Such ob- 
jects may be important high-frequency sources for LIGO 
and other ground based gravitational wave detectors. 
Any future detection of such inspiraling binaries may 
therefore carry information about the internal structure 
of the neutron star and its equation of state . 

Considerable effort has gone into the theoretical mod- 
eling of binaries containing either two black holes (e.g. 0, 
SU) or two neutron stars (e.g. see also the 

reviews HHH). For the most part, BHNS binaries have 
been studied only in Newtonian gravity (but see 0]). 
Several authors [lj, [l5| have modeled BHNS binaries as 
Newtonian ellipsoids around point masses, generalizing 
the classic Roche model for incompressible stars (see, 
e.g. 16]) to compressible configurations. These ellip- 
soidal calculations have also been generalized to include 
relativistic effects 

E3 El EH Equilibrium models of 
Newtonian BHNS binaries also have been constructed by 
solving the exact fluid equations numerically and again 
treating the black hole as a point mass [2(|. Dynami- 
cal simulations of Newtonian BHNS binaries, including 
coalescence, merger, and the tidal disruption of the neu- 
tron star, have been performed, for example, by Q, l2ll 
and references therein] . Dynamical simulations have also 
been performed for the tidal disruption of normal stars 
orbiting massive black holes (see, e.g. and references 
therein). But, at best, these simulations typically em- 
ploy a relativistic treatment of the hydrodynamics with 
an approximate treatment of the gravitational field of 
the star and/or the gravitational tidal field of the black 
hole companion. There is greater urgency in solving this 
problem more carefully, now that the possible detections 
by the Chandra X-ray Observatory of tidal disruptions 
of normal stars by supermassive black holes have been 
reported [23| • Clearly, accurate and reliable models of bi- 
naries containing black holes, including BHNS binaries, 
require a fully relativistic treatment. 

This paper is the first in a series on quasiequilibrium 
models of BHNS binaries in full general relativity. Our 
models serve a dual purpose: (1) constant rest-mass, 
quasiequilibrium sequences parametrized by binary sepa- 
ration approximate the adiabatic inspiral phase and can 
be used to locate the ISCO or the onset of tidal insta- 



bility, and (2) individual models provide numerically ex- 
act initial data for future relativistic dynamical simu- 
lations. Similar to earlier treatments of binary neutron 
stars (e.g. [EGHEEIEJ) an d binary black holes (e.g. 00, 
we adopt the conformal thin-sandwich decomposition |2 1| 
of the constraint equations of general relativity. We solve 
these equations for the gravitational field, together with 
the relativistic Euler equation for the matter, assuming 
the presence of a helical Killing vector. For the back- 
ground geometry, we adopt the metric corresponding to 
a single black hole in Kerr-Schild coordinates, thereby 
accounting for the black hole companion of the neutron 
star. 

In this first paper we provide solutions for the simplest 
case in which the mass ratio is large, i.e., Mbh 3> M^s- 
While the estimate above shows that, in this limit, the 
internal structure of the neutron star and tidal distortion 
is small, it allows us to formulate the BHNS problem in 
full general relativity and to construct numerical models 
of a neutron star in the presence of a black hole compan- 
ion. We also adopt a polytropic equation of state for the 
neutron star matter. In a future paper we will relax these 
approximations and will construct binaries with compan- 
ions of comparable mass, for which the internal struture, 
tidal distortion and tidal disruption are potentially very 
important. 

The paper is structured as follows. In Section ^ we 
describe how the BHNS problem can be solved in New- 
tonian gravity. A very similar approach is adopted in 
the relativistic treatment, so it is useful to outline this 
approach first in the much more transparent Newtonian 
framework. In Section IlIII we develop the relativistic for- 
malism. We introduce the assumption of extreme mass 
ratios in Section llVl In Section we present numerical 
results, including code tests and constant rest-mass se- 
quences. We provide a brief summary in Section IVII A 
detailed description of the numerical strategy that we 
adopted in solving the relativistic equations, together 
with Kerr-Schild background expressions that appear in 
these equations, can be found in Appendix lAl Through- 
out this paper we continue to adopt geometrized units 
(G = c=l). 



II. BHNS BINARIES IN NEWTONIAN 
GRAVITY 

Before turning to the relativistic problem, it is useful 
to lay out in this Section how BHNS binaries can be 
constructed in the much simpler Newtonian framework 
(compare |2(j). In the Sections that follow we will adopt 
a very similar strategy and algorithm for the construction 
of the same binaries in general relativity, and will refer 
back to this Section for clarity. 
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FIG. 1: Coordinate system that we use for the construction 
of BHNS binaries. The origin is centered on the neutron star, 
the black hole is located at the position ibh, and the axis of 
rotation, which is parallel to the z-axis, is located at x to t- For 
the extreme mass ratio Mbh S> Mns adopted in this paper 
we have ibh = x la t- 



A. The Poisson Equation 

In Newtonian gravity, the black hole can be described 
by a point mass Mbh which gives rise to the gravitational 
potential 



Mbh 
rBH 



(2.1) 



where tbh is the distance from the point mass. The total 
gravitational field <p can be written as a sum 



!>BH 



0NS, 



(2.2) 



where 0ns is the neutron star contribution to the poten- 
tial. Since 0bh is a homogeneous solution to the Poisson 
equation (for tbh > 0), we can write the Poisson equa- 
tion as 



D 2 (f> = D 2 (j) NS = 4np 
where D 2 is the Laplace operator. 



B. The Integrated Euler Equation 



(2.3) 



We assume a polytropic equation of state 
P 



Kp 



(2.4) 



where P is the pressure, po the (rest-mass) density, n 
the polytropic index and k the polytropic constant. For 
corotating equilibrium configurations the Euler equations 
can then be integrated analytically to yield 



P 



1 



( n + l) -+cj>- -fi 2 ((. 



Po 



) 2 +y 2 )=C, (2.5) 



where f2 is the orbital angular velocity, C is a constant 
and where we have assumed that the star orbits about an 
axis parallel to the z-axis located at x = x ro t and y = 0. 
In the following we will also assume that the black hole 
and neutron star are aligned with the x-axis (see Fig. ^| . 

We point out that it is unlikely that the neutron star 
viscosity could be strong enough to lock the star into 
corotation during the binary inspiral |25l I2(t| . There- 
fore it would be more realistic astrophysically to assume 
an irrotational rather than a corotational fluid flow (see, 
e -g- ; and references therein). However, the numeri- 
cal implementation of irrotational fluid flow is much more 
involved than that of corotation, and it is therefore rea- 
sonable to focus on corotating configurations first (as it 
was done for binary neutron stars). 

Equations (|2.3(l and (|2.5J) are the fundamental equa- 
tions describing a stationary star orbiting a black hole in 
circular co-rotation. The solution will depend on three 
parameters, which can be chosen, for example, to be the 
neutron star mass 



Mi 



NS 



p d x, 



(2.6) 



the black hole mass Mbh and the separation d. For com- 
putational reasons we will use the maximum density of 
the neutron star as one of the three free parameters in- 
stead of the neutron star mass. Given these parameters, 
a solution to the problem will also provide the eigenvalues 
il and C, as well as the size of the neutron star. 



C. Elimination of Dimensions 

It is convenient to introduce the dimensionless density 
parameter 



q = 



in terms of which we have 



P 
Po' 



Po 



K~ n q n 



and 



P = K- n q n+1 . 



(2.7) 



(2.8) 



(2.9) 



Since k™/ 2 has units of length, we can also introduce 
dimensionless coordinates x = n~ n / 2 x, and similarly for 
y and z. The Laplace operator then scales as D 2 — K n D 2 , 
and the angular velocity as Cl = K n / 2 f2. In terms of these, 
equations (|2.3(l and l|2.5() become 



D 



Anq r< 



(2.10) 



and 



(n + l)q + <f> - -n 2 ((x - x rot ) 2 + f) = C. (2.11) 
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The integrated Euler equation (|2.11|) can also be rewrit- 
ten as 

q = ^-jry (c - <f) + ((x - x mt ) 2 + f)^J 



71 + 1 



(C - $ cff ) , 



where the effective potential $ c ff combines the gravita- 
tional and centrifugal potentials, 



$off = <t> - \fl 2 ((a; - Xrot) 2 + y 2 



(2.13) 



This demonstrates that, up to a constant, the density 
parameter q is proportional to the effective potential. 
Rescaling M = KT n l 2 M we also find 



M NS = / q n d 3 x 



(2.14) 



D. Rescaling 

As we will see below, it is convenient for numerical 
purposes to have the surface of the neutron star intersect 
the x-axis at fixed coordinate locations, say xa = — 1 
and xb = 1- This means that we have to rescale the co- 
ordinates with respect to the physical, non-dimensional 
size of the neutron star, say f e . For spherical stars, f e is 
the radius of the star, and otherwise it is half the star's 
diameter along the x-axis. We denote these new coor- 
dinates with hatted quantities, i — x/f e , and similarly 
for y and z. The derivative operator scales according to 
L>i = f e Di, and we also have = f e f2. The Poisson 
equation l|2.10|l then becomes 



D 2 NS = Anf 2 e q n 



(2.15) 



and the integrated Euler equation l|2.11|l 



(n + l)q + 4>- -Cl 2 ((x - i rot ) 2 + f) = C. (2.16) 

In these coordinates, the physical, dimensionless neutron 
star mass Mns can be computed from 



M NS = 



% d x = 



l d A x. 



(2.17) 



From the Poisson equation (|2.15l) it is evident that the 
gravitational potential <f> will depend on f e . To make 
this dependence explicit, we could introduce a rescaled 
potential 

(hs=rlh*s- (2.18) 

This expression can then be inserted into (j2.16|l . explic- 
itly introducing f e into the integrated Euler equation. 
We have yet to determine the scaling of the black hole 
contribution </>bh 
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(2.19) 



Here Mbh is the dimensionless black hole mass that we 
want to keep fixed during the iteration. 

With these choices, the integrated Euler equation 
(12.1611 can now be rewritten as 



(2.12) (n + l)q + (j)_ 



M BH 



NS 



l -Cl 2 ((X 



Xrot ) 



y 



c. 



^BH 

(2.20) 

This equation has to be solved together with the Poisson 
equation i|2.15[) . A self-consistent solution can be con- 
structed iteratively, as outlined immediately below. 



E. Iteration Scheme 

Instead of fixing the neutron star mass, we fix the max- 
imum density on the x-axis q max . A neutron star of a de- 
sired mass can later be constructed iteratively by varying 
<Zmax- We also fix the separation xbh in terms of the neu- 
tron star size r e , £bh = ^BH/^e- In the following we will 
locate the point mass at (— £bh,0, 0). The three input 
parameters are therefore q maxi £bh and Mbh- 

In this paper we also assume that the mass of the black 
hole is much greater than that of the neutron star 



Mbh > Mns > 



(2.21) 



in which case the rotation axis will go through the center 
of the black hole 



Zrot = 2:bh- 



(2.22) 



This assumption leads to an error in the location of the 
rotation axis of the order of Mns /Mbh times the binary 
separation; this becomes negligible as the mass ratio de- 
creases. Without this assumption, the location of the 
rotation axis i ro t has to be found as part as the iteration 
scheme. Various approaches can be chosen to identify 
x ro t; in Newtonian theory it coincides with the center of 
mass. 

In the limit of extreme mass ratios we expect that the 
neutron star will affect the physical solution only in a 
neighborhood of the star itself, so that we can confine the 
numerical grid to a finite domain around the neutron star 
that avoids the black hole. This significantly simplifies 
the problem, especially in the relativistic case, where the 
black hole interior does not have to be excised from the 
numerical grid. 

The iteration scheme starts with an initial guess for the 
density q, confined between xa and xb, and the star's 
physical size f e . 



1. Solution of the Poisson equation 

Given a density distribution q we can solve the Poisson 
equation l|2.15|l to find the neutron star potential </>ns- 
The elliptic equation is solved with PETSc algorithms 
that are described in more detail in the relativistic con- 
text below. 
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2. Determination of Eigenvalues 

Before the integrated Euler equation (|2.2U|) can be 
solved for the new matter distribution q, the eigenval- 
ues Cl, f e and C have to be determined. This can be 
done by evaluating the integrated Euler equation at the 
three points Xa, %b and xc- The latter is defined as the 
point along the x-axis at which q takes its maximum. At 
that point, we set q = g ma x, and at the other two points, 
which lie on the surface, we have q — 0. The potential 
</>ns is interpolated to these three points, yielding <j>^ s , 
4>§ s and </>ns- As we have discussed above, the neutron 
star potential <^>ns depends implicitly on f e . To make 
this dependence explicit, we compute 0ns = ^Ns/^e at 
the three points, which finally results in the following 
three equations 

° — r e ™S - -A o 12 \ XA Xrot ) 
r e r BH Z 

C = f 2 J* s - - \&{x B - x mt ) 2 (2.23) 

C = fg^NS ~ - '^c _ \^ 2 ( X C - ^rot) 2 + (n + l)«max- 

These equations can be solved iteratively for the eigen- 
values O, f e and C. Without the assumption of extreme 
mass ratios a fourth condition has to be added so that 
x ro t can be determined together with the three other 
eigenvalues. 

3. Rescaling of neutron star potential 

With the new value of f e , the rescaled neutron star 
potential is computed from 

(f>NS=f 2 J NS - (2.24) 

4- Solution of the integrated Euler equation 

Given the three constants (l, f e and C the integrated 
Euler equation H2.20fl can now be solved everywhere. The 
innermost closed surface centered on the origin on which 
the density q equals zero is identified as the stellar sur- 
face. Immediately outside this surface, q derived from 
equation (|2.2U|) becomes negative; once we identify the 
surface, we properly set q to zero everywhere outside. 
After obtaining the solution of the integrated Euler equa- 
tion, we evaluate the residual of the Poisson equation. If 
this residual is smaller than a specified tolerance, the it- 
eration is terminated; otherwise another full iteration is 
performed. 

This completes the construction of a Newtonian BHNS 
binary for a given black hole mass Mbh, a binary sep- 
aration £bh and a maximum neutron star density q max . 
Binaries with neutron stars of a given mass can then 



be constructed iteratively by varying g max , and constant 
mass sequences can be computed by varying xbh- 



III. BHNS BINARIES IN GENERAL 
RELATIVITY 

Constructing fully relativstic BHNS binaries amounts 
to finding a spacetime metric 

ds 2 = g ab dx a dx b = -a 2 dt 2 +~ ftJ (dx 1 + /3 i dt)(dx j + ft dt) 

that solves Einstein's equations together with a matter 
distribution 

Tab = {po + Pi + P)u a u b + Pg a b (3.2) 

that satisfies the relativistic equations of hydrodynam- 
ics. In the above "ADM" (2?| form of the metric, a is 
the lapse, (i % is the shift vector, and 7y is the spatial met- 
ric. In the stress-energy tensor, pi is the internal energy 
density and u a the fluid four-velocity. 

In this paper we focus on the construction of quasiequi- 
librium initial data, so that we only seek to find a solu- 
tion to the constraint equations of Einsteins equations. 
We solve these constraints within the conformal thin- 
sandwich decomposition, which provides a natural frame- 
work for the construction of quasiequilibrium (Section 
MI Al compare [Til II2I l24jl. As in Newtonian gravita- 
tion, the assumption of equilibrium allows us to model 
the star as a solution to an integrated Euler equation 
( Section IIII Bp . The numerical implementation of these 
equations again requires a rescaling, which we discuss in 
Sections (fTTTT^f) and PTTTj)l . 

The constraint equations fix only some of the gravi- 
tational fields, while others, which determine the "back- 
ground" geometry, can be chosen freely, but according 
to the astrophysical context. In our approach to con- 
structing BHNS binaries, we choose for the background a 
Schwarzschild black hole in Kerr-Schild coordinates (Sec- 
tion llll E|l . Our background therefore describes the black 
hole in the binary and plays a role that is equivalent to 
the analytic black hole potential 0bh in the Newtonian 
context. We then solve the constraints together with the 
integrated Euler equation to model a neutron star in the 
presence of this black hole. 



A. The conformal thin-sandwich equations 

The conformal thin-sandwich decomposition has been 
used extensively for the construction of binary neutron 
stars as well as binary black hole systems (e.g . |5|,|6ll3)' 
and has been developed independently in |2J|. We refer 
to the reviews ^3 ^2 f° r a detailed derivation, and only 
state the most important results here. 

The spatial metric 7^ is conformally decomposed into 
a conformally related background metric 7^ and a con- 
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formal factor i[>, 



7ij = ■ 



(3.3) 



Below we will assume that 7y is given by a Kerr-Schild 
metric. For the construction of quasiequilibrium data in 
a corotating frame it is natural to assume that the time 
derivative of the conformally related metric vanishes, 



u ij = dtlij = 0. 



(3.4) 



For the trace of the extrinsic curvature we choose that of 
the Kerr-Schild background, and set its time derivative 
to zero. With this choice all the freely specifiable quan- 
tities are fixed, and Einstein's equations can be used to 
determine the lapse a, the shift /3 1 , and the conformal 
factor tp. 

The trace-free part of the evolution equation for the 
spatial metric yields 



lb 6 - 



(3.5) 



where A lJ is the trace-free part of the extrinsic curvature 
, — ip w A i: > is the conformally related trace-free 
extrinsic curvature, and where we have used Uij = 0. The 
longitudinal operator L is defined as 

(L(3Y j = D l (3 3 + & j3 i - \f j D k (3 k . (3.6) 

o 

Equation (|3.5|l can now be inserted into the momentum 
constraint, which yields an equation for the shift (3 l 



A L p l -(Lp) l W J \n(aij- 6 )--aD i K = lQita^f. (3.7) 

o 

Here j a = —"f a b n c T bc is the momentum density observed 
by a normal observer and the vector Laplacian is 

A L (3 l = Dj{Lpji = D 2 (3 l + \&bj^ + R\(3\ (3.8) 

3 J 

where R l j is the Ricci tensor associated with the confor- 
mally related metric 7^. 

The conformal factor tp is then determined from the 
Hamiltonian constraint 

b^- l^R- ^ 5 K 2 + l^Aij Av=-2^p, (3.9) 

where p = n a ribT ab is the energy density observed by a 
normal observer and where 

b 2 ip = fibib^ = - fV,i (3.io) 

is the Laplace operator associated with 7^ . 
To derive an equation for the lapse we assume 

d t K = (3.11) 



(but note that Using the trace of the evolution 

equation for the extrinsic curvature we find that this con- 
dition implies 



D 2 a = a{K i0 K i] + 4ir(p + S)) + l3 l DiK 



(3.12) 



= aff'M,,^ + -K 1 + Att( P + S)) + I3 1 D,K, 

where S = 7 T a b. This equation involves the Laplace 
operator with respect to the physical metric 7ij , which we 
do not know a priori. It is therefore convenient to rewrite 
this operator in terms of the conformal background met- 
ric 7y. To do so, we combine equation l|3.12|l with the 
Hamiltonian constraint (|3.9() . which yields 



1 -aip-'' 'A l0 A 11 + ^-aip 5 K 2 + \aipR 
12 8 



+^ 5 p l b t K + 2naij 5 (p + 2S). (3.13) 



Equation (|3.13|) is the generalization of the dtK = con- 
dition for a non-zero (but constant in time) K. It reduces 
to the more familiar maximal slicing equation for K = 0. 

Equations (|3.9|) , (|3.7|) and l|3.13|l now form the confor- 
mal thin-sandwich equations for the unknown gravita- 
tional field quantities tp, (3 l and a. These equations are 
equivalent to Poisson's equation (|2.3|l in the Newtonian 
framework. 



B. The integrated Euler equation 

For corotating equilibrium solutions, the rclativistic 
Euler equation can be integrated analytically to yield 



= 1 + C, 



(3.14) 



where C is a constant (we add one to the right hand 
side so that C reduces to the equivalent constant in the 
integrated Newtonian Euler equation (|2.5|) ~). and where 
h is the specific enthalpy 



h — exp 



dP 



Po 



P 



(3.15) 



(see, e.g., yjj.) Here pi is the internal energy density. 
Note that the integrated Euler equation is consistent with 
von Zeipel's theorem, which states that in uniformly ro- 
tating, perfect fluid stars, surfaces of contant P, po and 
u* all coinide. 

The time component of the four-velocity u* can be ex- 
pressed in terms of the relative velocity v between the 
matter and a normal observer as 



1 



( 1 _ u 2)l/2- 



(3.16) 



In a corotating coordinate system, where u % ~ 0, u* can 
be found from the normalization condition u a u a = — 1, 
which yields 



u* = (a 2 -lijPP) 



-1/2 



(3.17) 



7 



This expression, together with the enthalpy (|3.15() , can be 
inserted into the integrated Euler equation l|3.14[) . yield- 
ing a relation between the fluid and gravitational field 
variables. 

As before, we adopt a polytropic equation of state (|2.9|1 
and introduce the density parameter q (|2.7|1 . The en- 
thalpy (|3.15|) can then be written as 

h= Po+Pi + P =l + (l + n)q, (3.18) 
Po 

where we have used 

Pl =nP = nn- n q n+1 . (3.19) 
The integrated Euler equation then becomes 

(3.20) 



1 



(„*(! + O-l) 



In the Newtonian limit this equation reduces to the inte- 
grated Newtonian Euler equation (|2.12|) . 



D. Rescaling 

As in the Newtonian problem, we would like the star's 
surface to intersect the x-axis at constant coordinate val- 
ues of xa = — 1 an d £g = 1 . This again requires a rescal- 
ing with respect to the star's dimension, which we again 
call f e (see Section Hi Dll . Denoting these rescaled quan- 
tities with hats, e.g. D = f e D, R = f&R, A*3 = r e A ij , 
K = f e K, we find the Hamiltonian constraint 



1 , a . 1 
12 



1 



D 2 ^ = -i{jR+—ilj 5 K 2 --ij- 7 A i3 A tj -2TnP 5 f 2 e p, (3.23) 



the Momentum constraint 



k L ft = (LffiDj ln(c«/r 6 ) + -aD l K + lGwa^fl?, 

o 



and the lapse equation 



(3.24) 



8 12 8 

+^ 5 f3 t D t K + 2Trai; 5 f 2 e {p + 2S). (3.25) 



C. Elimination of Dimensions 

As in Section (|II C|l for the Newtonian treatment, we 
can now eliminate dimensions by rescaling all dimen- 
sional quantities with respect to appropriate powers of 
k™/ 2 , which has units of length. In particular we have 
D = K n / 2 D, R = K n R, K = K n / 2 K, and A lj = n n / 2 A l K 
Denoting these dimensionless quantites with bars, we ar- 
rive at the following system of field equations 

ib 6 - ■■ 
A*> = ^-(L/?) y 
2a 



D 2 ^ = -^R + —^K 2 - ^- 7 A l3 A lj - 2n4> 5 p 



A L f3 l = (L(3) ij D j ln(ai(j- 6 ) + -aD l K +16™^? 



D 2 (aiP) = -a^- 7 A i:j A ij + \a^p b K 2 + laipR 



lJJ - ' 12' 

1.5 ai n. i> i o„„,„/,5 / 



+ip b i3 l D l K + 2 7 ra#'(p + 2S). 



(3.21) 



Here the matter sources can be expressed in terms of the 
dimensionless quantity q 



.., ,1 + (1 + n)q 
p = K n p = q n ( / 2 - q 



J = K'*3 = q 



1-v 2 
„ 1 + (1 + n)q 
1 — v 2 a 



(3.22) 



§ = K " S = g» ( 1 + a + n) % 2 + 3q 



The field equations l|3.21|l have to be solved together with 
the integrated Euler equation l|3.20|l . 



In the above equations, the extrinsic curvature is com- 
puted from 



ib 6 . .. 
A« = I-(L/3) y . 
2a 



(3.26) 



Inserting H3.17fl into l|3.20|l we also find the integrated 
Euler equation 



l + C 



1 + n \ {a 2 -iWypp) 1 / 2 



(3.27) 



We note that the metric is dimensionless, so that 7^ = 
lij = lij- For ease of notation, we decorate the confor- 
mally related metric with a hat. 



E. Kerr-Schild Background 

As discussed above, we account for the presence of the 
black hole by choosing a background metric that, in the 
absence of the neutron star, describes a Schwarzschild 
black hole. This black hole can be expressed in differ- 
ent coordinate systems. Probably the simplest choice 
would be standard isotropic coordinates (in which case 
the presence of the black hole is actually absorbed in the 
conformal factor, while the background 3-metric is flat). 
The disadvantage of these coordinates is that they do 
not penetrate into the black hole interior. We expect that 
such penetration is crucial when dealing with companions 
of comparable mass, in which case the spacetime in the 
vicinity of the black hole and its apparent horizon have to 
be solved for. Painlcvc-Gullstrand coordinates are both 
remarkable simple and regular on the black hole horizon 
(see [2!| for a recent discussion). However, these coor- 
dinates lead to spatial slices that are not asymptotically 
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flat, complicating the definition of mass and angular mo- 
mentum. We have therefore chosen to adopt Kerr-Schild 
(or ingoing Eddington-Finkelstein coordinates), in terms 
of which the Schwarzschild metric for a single, nonrotat- 
ing black hole can be written as 

ds 2 = g ab dx a dx b = (r] ab + 2Hl a l b )dx a dx b (3.28) 

(compare HHUEilll.) Here 

^ = MBH = iWBH = iMBH ) 

t"bh r B H r e r B H 

is the analogue of the Newtonian potential l|2.19[l and 
scales in exactly the same way. As in the Newtonian 
case 

1 /2 

rBH = [(x - x B n) 2 + (y- ma) 2 + (z— z B h) 2 ] 

(3.30) 

is the coordinate distance to the center of the black hole 
at the coordinate location (— xbh,0,0) (rBH is the famil- 
iar areal radial coordinate centered on the Schwarzschild 
black hole). The vector l a is null with respect to both 
g a b and the "Minkowski" metric r\ ab and can be written 
as 



1. 



h = V 



X" 

r-BR ' 



(3.31) 



Comparing the Kerr-Schild metric (|3.28|l with the 
ADM form of the metric H3.1(l shows that we can identify 
the lapse as 



a w = (1 + 2H) 



-1/2 



the shift as 



and the spatial metric as 

7y = ihj + '2111,1 r 



(3.32) 
(3.33) 
(3.34) 



Here we use the subscript BH to distinguish these back- 
ground quantities from the lapse a and shift (3 % of the full 
solution. 

We note that l a is a four-vector, which is therefore 
raised with the spacetime metric 



l a = (r) ab - 2Hl a l b )l b = r) ab l b 



(3.35) 



as opposed to the spatial metric. The inverse of the spa- 
tial metric is therefore given by 



f j = «bh ((1 + 2H W 3 ~ 2HlH j ) 



(3.36) 



Other useful quantities that appear in many of the equa- 
tions are 



T' 
K 



BH 



rBH 



■(2P m -m 3 i k ) 



(3.37) 



tbh 
rBH 



y k £) k = "22^(3 + 8H)P (3.38) 



(3.39) 



Several other quantities will be required for the solution 
of the thin-sandwich equations in a Kerr-Schild back- 
ground, and we will provide those as they appear in the 
equations. 

Just like the Newtonian black hole potential i|2.1[l is a 
vaccum solution to the Poisson equation (|2.3J) . the back- 
ground lapse (|3.32|) . shift (|3.33|) and the trivial back- 
ground conformal factor V'bh — 1 solve the field equa- 
tions (|3.23(l - H3.25|) identically in the absence of mat- 
ter sources, p — f = S = 0. In further analogy with 
the Newtonian problem we therefore solve the relativis- 
tic equations by computing corrections to the background 
quantities that account for the presence of the neutron 
star. Details of the numerical strategy and implementa- 
tion are presented in Appendix 1X1 

Once a model has been constructed, we can compute 
the neutron star's rest mass from 



Mo = J pou a d 3 Y, a — J pou t \/—gd 3 x. 



(3.40) 



In terms of the variables used in our code this can be 
rewritten as 

M = «T"/ 2 M = fl { aip 6 u t y /j B -^q n d 3 x, (3.41) 



where 7bh is the determinant of the spatial Kerr-Schild 
background metric. 



IV. EXTREME MASS RATIOS 

In this paper we focus on the construction of neutron 
stars in the presence of a black hole companion by as- 
suming extreme mass ratios Mbh >> Mns- This ap- 
proximation simplifies the problem in a number of ways. 

Most importantly, we may assume that the spacetime 
is affected by the presence of the neutron star only in a 
neighborhood of the neutron star. The numerical grid 
therefore needs to cover only a region around the neu- 
tron star, and does not need to encompass the black hole. 
That means that we do not need to excise the black hole 
interior, simplifying the computational domain and elim- 
inating the need for interior black hole boundary condi- 
tions. This also means that the computational grid can 
be chosen significantly smaller, which speeds up numeri- 
cal calculations. 

Another simplification arises from the fact that in the 
extreme-mass-ratio limit, the binary's center of mass re- 
sides at the center of the black hole, £ ro t = £bh, which 
eliminates the need to locate i ro t iteratively. In the gen- 
eral case, the radius a; ro t can be identified, for example, 
by noting that the solution becomes quasistationary in 
the corotating frame and therefore must satisfy the addi- 
tional condition that Madm = Mpcnmar when the center 
of mass is correctly located (see |33| for an illustra- 
tion). Finally, for the construction of constant-mass se- 
quences, the irreducible mass of the black hole has to 
be held constant, which requires an additional iteration. 



9 




0.1 0.2 0.3 0.4 

q 



FIG. 2: Rest mass Mo of n = 1 polytropes in isolation as 
function of the density parameter q. The crosses denote the 
exact results for spherical, non-rotating neutron stars, and 
the lines correspond to numerical results for a fine grid (16 
gridpoints across a neutron star radius; denoted by the solid 
line), a medium grid (8 points; dashed line) and a coarse grid 
(4 points; dotted line). In the inset we show the rescaled 
numerical error, i.e. the difference between the fine-grid solu- 
tion and the exact solution, multiplied by a factor of 16 and 
the difference between the medium grid solution and exact 
solution, multiplied by factor of 4. The convergence of the 
rescaled error functions establishes second-order convergence 
of our implementation. 

This iteration is also unnecessary in the extreme mass 
ratio limit. 



V. NUMERICAL RESULTS 

In this Section we present numerical results, both for 
neutron stars in isolation (as a code check, Section IV A|) 
and for binary sequences of variable separation but con- 
stant rest mass (Section IV til . We discuss tidal disrup- 
tion in Section IV CI For all results in this Section wc 
adopt a polytropic index of n = 1. 



A. Neutron Stars in Isolation 

As a simple code test we can use our formalism to con- 
struct neutron stars in isolation (Mbh = 0), for which the 
solution is given by the non-rotating, spherical Tolman- 
Oppenheimer-Volkoff (TOV) solutions TOV solu- 

tions can be considered "exact" since they satisfy a set 
of ordinary differential equations that essentially can be 
solved to arbitrary accuracy. 
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FIG. 3: Same as Fig. [5] except for the isotropic neutron star 
radius R. The three different resolutions are almost indistin- 
guishable and agree well with the exact solution marked by 
the crosses. However, the inset with the rescaled numerical er- 
rors shows that the coarse resolution with only 4 grid points 
across the neutron star radius is not yet in the convergent 
regime. 

In Figs. and we show the exact values of the rest 
mass Mq and the isotropic radius R as a function of the 
density parameter q, together with our numerical results 
for three different grid resolutions. For these tests we 
imposed the outer boundaries at X out = 2 (see Section 
IA l|l . i.e. at twice the neutron star radius (see Sections 
III PI and IIII DD . We adopted numerical grids of (64 x 
64 x 32), (32 x 32 x 16) and (16 x 16 x 8), corresponding 
to 16, 8 and 4 gridpoints along a neutron star radius. 
Focussing on the rest mass Mq plot, it can be seen how 
the numerical results converge to the exact results. In the 
insets we also show the rescaled numerical errors, which 
establish second order convergence of our code. We find 
very similar results for the gravitational (ADM) mass M. 

B. Constant Rest Mass Sequences 

Inspiraling BHNS binaries conserve the rest mass of 
the neutron star and the irreducible mass of the black 
hole. We construct sequences of constant neutron star 
rest mass by iterating, for a given binary separation xbh, 
over the maximum density parameter g max until a neu- 
tron star of the desired rest mass has been found to within 
a specified tolerance. 

Keeping the black hole's irreducible mass constant 
would involve evaluating the area of its event horizon (or, 
in practice, its apparent horizon). The difference between 
the black hole background mass A/bh and the black hole's 
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FIG. 4: The orbital frequency as a function of binary separa- 
tion for mass ratios Mbh /Mns = 100 and three different neu- 
tron star rest masses Mns = 0.05, 0.1 and 0.15. Our results 
typically agree with the Kepler frequency 15. H to within one 
percent. The dot marks the Kepler frequency D.A4bh = 6 _3//2 
at the innermost stable circular orbit at ibh = 6A/bh. For 
these simulations we used a grid of (64 x 64 x 32) gridpoints 
and imposed outer boundaries at X ut = 4. 



irreducible mass is in the order of the binary's binding 
energy [HI . In this paper we neglect effects of the neu- 
tron star on the black hole, and hence neglect this change 
in the black hole mass. Strictly speaking, we therefore 
construct sequences of constant neutron star rest mass 
Mo and black hole background mass A/bh- 

In the test-mass limit, the orbital frequency f£ is given 
by the relativistic version Kepler's third law 



BH 



-3/2 



(5.1) 



FIG. 5: Contours of the density parameter q (solid lines) and 
the right hand side of the integrated Euler equation 13.2UH 
(dotted lines) for Mns = Mbh/10 = 0.05 at a separation 
of ibh = 12.6A4bh (the center of the black hole is located at 
x — —5.0 and y = 0). The innermost contour has a density ra- 
tio q/f/max = 10/11, and successive contours have q decreasing 
by the same decrements Ag/q max = 1/11. The saddle point 
between the neutron star and the black hole defines the inner 
Lagrange point L\ and is marked by the cross. The equipo- 
tential surface passing through L\ is the relativistic Roche 
lobe. The proximity of L\ to the neutron star surface indi- 
cates that the neutron star is close to being tidally disrupted. 
The saddle point "above" the neutron star is the outer La- 
grange point L<2- The inset shows both q (solid line) and the 
right hand side of 1)3. 20^ (dotted line) along the x-axis as a 
function of x. The saddle point is close to the minimum of 
the dotted line at about x = —1.3. 
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In Fig.0]we plot J7A/bh for mass ratios A/bh/A/ns = 100 
and three different neutron star rest masses Af NS =0.05, 
0.1 and 0.15 (the maximum allowed rest mass for isolated, 
n = 1 polytropes in these units is A/ns = 0.18.) The 
corresponding compaction M^sl ' Raxcai of these neutron 
stars in isolation are 0.042, 0.088 and 0.145 (here i? a rcai is 
the neutron star's areal radius). All frequencies typically 
agree with the relativistic Kepler frequency to within bet- 
ter than one percent, and, as expected, the frequency is 
independent of the neutron star mass. The dot in Fig. 0] 
marks the orbital frequency OA/eh = 6~ 3 / 2 = 6.804 at 
the innermost stable circular orbit at Xbh = 6 A/bh • 



TABLE I: Numerical values for a constant rest-mass sequence 
with Mbh = 0.5 and Mns = 0.05 close to the onset of tidal 
disruption at a separation of ibh = 11.9Mbh (see Fig- - 



C. Tidal Disruption 

In corotating binaries, tidal disruption occurs when 
the neutron star overflows its Roche lobe. In Newto- 
nian physics, the Roche lobe is defined as the innermost 
equipotential surface of the effective potential (|2.13|l that 
encompasses both binary companions. In a contour plot 
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FIG. 6: Same as Fig. [5] except for a separation of Sbh = 
11.9Mbh (the center of the black hole is located at x — —4.18 
and y = 0). The Lagrange point L\ is now located on the 
star's surface, indicating the onset of tidal disruption. 



in the equatorial plane the Roche lobe forms a "figure 
eight" , in which each teardrop-shaped region contains 
one companion. The equipotenial surface pinches off at 
the saddle-point of the potential, which is called the in- 
ner Lagrange point L\. When the neutron star overflows 
its Roche lobe, matter flows across this inner Lagrange 
point and accrets onto the black hole. 

By virtue of the integrated Euler equation (|2.12(l the 
density parameter q is proportional, up to a constant, to 
the effective potential $ c ff ■ The Roche lobe can therefore 
be constructed equivalently from contours of the right 
hand side of the integrated Euler equation (|2.12l) . This 
immediately suggests the construction of a relativistic 
Roche lobe from the right hand side of the integrated rel- 
ativistic Euler equation (|3.20|1 . For a relativistic equilib- 
rium configuration, the Roche lobe and Lagrange points 
can therefore be identified in complete analogy to the 
Newtonian case. 

In Figs. and 10 we show contours for two config- 
urations with Mns = Mbh/10 = 0.05 in the equatorial 
plane. Table [i] also provides numerical values for this 
sequence. Fig. (0 shows the binary at a separation of 
xbh = 12.6Mbh (the second line in Table QJ. For this 
calculation and all others in this Section we used a grid 
of (48 x 48 x 24) gridpoints and imposed outer boundaries 
at X out — 3. Density contours are marked by solid lines, 
while contours of the right hand side of (|3.20|) . i.e. con- 
tours of the relativistic effective potential, are marked by 
the dotted lines. The inner Lagrange point L\ is marked 
by the cross. The contour passing through this point 
is the Roche lobe. The outermost solid contour, which 



FIG. 7: The tidal separation Stid/^NS as a function of mass 
ratio Mns/Mbh- The crosses mark our numerical results for 
A'/ns = 0.01, the dot marks a result obtained by Uryu & 
Eriguchi [gjj, and the solid line denotes the expected qual- 
itative scaling relation 11.311 . (Mns/A^bh) -1 ^ 3 - The dashed 
line, 2.35(Mns/Mbh) _1,/3 , provides an approximate fit to our 
numerical data. 
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38.5 


6.43 


0.00405 


0.1 


61.7 


5.18 


0.00199 


0.1 




5.17 


0.00208 



TABLE II: Numerical values for the onset of tidal disruption 
of an Mns = 0.01 neutron star (whose isotropic radius in 
isolation is .Rns = 1.23). The bottom line are the results 
from the Newtonian models of Uryu & Eriguchi [2fj| (see their 
Table 2). 



marks the surface of the star, is still inside the Roche 
lobe, so tidal disruption has not yet set in. 

In Fig. © we show the same binary, but at a slightly 
smaller separation of xbh = 11.9-Mbh (the last line in 
Table Pi. Now the Lagrange point L\, again marked by 
a cross, is located on the surface of the star, indicating 
that the star now fills out its Roche lobe, and that tidal 
disruption sets in. This illustrates how the onset of tidal 
disruption can be identified by the Lagrange point L\ 
reaching the surface of the star, and simultaneously the 
star filling out its Roche lobe. 

It is now of interest to see whether we recover the 
scaling (|1.3fl . To better compare with Uryu & Eriguchi 
[2fjj , who have located the tidal separation in Newtonian 
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BHNS binaries, we compute their measure of separation 

Stid = x B K + -rj— [ xpd 3 x (5.2) 
Mm J 

(denoted d G _in HJ.) In Table GU and Fig. □ we show 
results for M NS = 0.01, for which R NS /M NS = 123. 
For these configurations relativistic effects are sufficiently 
weak to make a comparison with Newtonian results 
meaningful. According to the scaling relation 1)1. 4[l . tidal 
disruption will occur over a reasonable range of mass ra- 
tios Mns/AIbh -C 1. Crosses in Fig. J7J) denote results 
from our code, and the dot marks the only result from 
H3 for n = 1 and A/ ns /M B h < 0.1 (see also the bot- 
tom line in Table ITl) . Our result agrees with that of [2(j 
to within a few percent (and even better for the tidal 
separation Stid), which is well within the error that one 
expects for a mass ratio of Mns/JWbh ~ 0.1 within our 
approximation Mns/A/bh *C 1. For this configuration 
both M NS /R N s ~ 0.01 and A/ B H/s t id « 0.01 so that rel- 
ativistic effects could also introduce a deviation of a few 
percent. 

The scaling of the tidal separation with the mass ratio 
is very close to that predicted by (|1.3(l : in fact, the crude 
estimate for the tidal separation differs by only about a 
factor of 2.35 from our numerical results. 



VI. SUMMARY 

We have developed a new relativistic formalism for the 
construction of BHNS binaries in quasiequilibrium. We 
solve the constraint equations of general relativity, de- 
composed in the thin-sandwich formalism, together with 
the integrated Euler equation for the neutron star mat- 
ter. The background geometry describes a Schwarzschild 
black hole, which accounts for the neutron star's binary 
companion. Our approach yields self-consistent solutions 
describing a neutron star and its self-gravity in the pres- 
ence of a black hole companion. 

In this paper we assume that the black hole mass is 
much larger than that of the neutron star, which sim- 
plifies the problem in a number of ways. In this case, 
the axis of rotation, which in general is located at the 
binary center of mass, passes through the center of the 
black hole. Also in this limit, we can assume that the 
neutron star only affects the spacetime in a neighbor- 
hood around the star, so that the numerical grid can be 
restricted to such a neighborhood and does not need to 
cover the black hole. This eliminates the need for black 
hole excision. Furthermore, it is appropriate to assume 
that the black hole's irreducible mass, which is essentially 
constant during a quasiequilibrium inspiral, is equal to 
the background mass, which eliminates the need for an 
iteration. 

Adopting this approximation, we have tested our code 
for isolated neutron stars and for constant rest mass se- 
quences, finding excellent agreement with Kepler's third 
law for the orbital frequency. We have also developed a 



fully relativistic framework for the identification of the 
Roche lobe and Lagrange points, which can be used to 
locate the binary separation at which tidal break-up sets 
in. For weakly relativistic n — 1 polytropes, for which 
the tidal separation is outside of the innermost stable cir- 
cular orbit even for large mass ratios, we find that tidal 
disruption occurs at about twice the separation suggested 
by the crude scaling relation Stid/^NS ~ (Mbh/Mns) 1 / 3 . 

Astrophysically more interesting scenarios require 
Mbh ~ Mns, for which the above assumption of ex- 
treme mass ratios has to be relaxed. We plan to provide 
solutions to black hole - neutron star binaries with com- 
panions of comparable mass in the near future. 
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APPENDIX A: NUMERICAL STRATEGY 

In this Appendix we provide details on the numerical 
strategy and implementation for the solution of the field 
equations l|3.23|l - (|3.25|) . together with the integrated 
Euler equation (|3.27|l . 



1. Numerical Grid 

The only symmetry assumption that can be made to 
simplify the problem is equatorial symmetry across the 
z = plane. In the Newtonian case, or in the case of 
maximal slicing K — 0, the solutions are also symmetric 
across the y — plane. But the presence of a non-zero K 
as in the case considered here there is no such symmetry 
condition across the y = plane (see also Appendix E 
in |36|). The numerical grid can therefore be chosen to 
cover a region x min < x < x max , y min < y < j/ max and 
< z < z max . For the results in this paper we chose 
the limits to be the same in all dimensions, — x m i n — 
x max = -ymin = Z/max = z max = A out . In this paper we 
also assume that Mbh 3> -Mns and restrict the numerical 
grid to a neighborhood of the neutron star, which requires 
A ou t < i£bh (the origin of the coordinate system is at the 
center of the neutron star, compare Fig. 
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2. Elliptic Solver 



where we have used the abbreviations 



A variety of different approaches can be chosen to 
solve the elliptic equations (|3.23|) - (|3.25|) . together with 
the algebraic equation l|3.27fl . We adopt finite differ- 
ence methods, and solve the resulting finite difference 
equations on uniform grids in a parallel environment us- 
ing DAGH software [3jj. Having adopted a Kerr-Schild 
background, the operators in equations l|3.23[) - (|3.25|) 
are not flat, which poses an additional challenge. One 
possible method of solving such a problem is to replace 
the "covariant Laplace" operator with a flat Laplace op- 
erator and solve the equation iteratively. In this paper 
we instead use PETSc software |38| to solve linear elliptic 
equations of the form 



* = 1 + * 



D 2 f + u f f = r f 
directly. Here D 2 is a covariant Laplace operator 



(Al) 



D 2 f = f j f,i j -t i f, i (A2) 

and Uf and r/ are functions. In the following we assume 
that the metric 7 y and the connection coefficients T 1 
are given by the background quantities (|3.36[) and 1)3. 38[) . 
Equations 1)3. 23|) - (|3.25|) can then be solved by casting 
them into the form (|A1|) . 



3. Hamiltonian Constraint 

We write the conformal factor as a sum of contributions 
from the black hole and the neutron star 



1p = 4>BH + IpNS- 



(A3) 



In the absence of the neutron star the background metric 
jij is already a solution to the constraints, so 



'0BH = 1 



(A4) 



identically. The split 1A3(I is convenient since i/'ns now 
satisfies a Robin fall-off conditions at large separation 
from the neutron star. 

Since the Hamiltonian constraint (|3.23|l contains non- 
linear terms, the equation has to be linearized before it 
can be solved with the linear operator IjAlfl . A solution 
to the non-linear problem is then constructed iteratively. 
We therefore write the solution ip^s as a sum of the pre- 
vious result ^ns an d a correction JV'nSj 



(A5) 



Inserting this into (|A3|) and l|3.23|l then yields, to first 
order in J^nS) 



(A6) 



= -D 2 ^ + -VR + —y 5 K 2 - -^- J A 2 - 2n^ 5 f 2 p 
8 12 8 



and 



NS 



A 2 — A- ■ AV 



(A7) 



(A8) 



Defining the residual of equation 1)3.23)1 as 



Res^ = D 2 ip - \ijiR - -^-ij) b K 2 + \i>-~' 'A 2 + 2-K^f\p 
8 12 8 



we can rewrite equation i)A6)l in the form i)Al)l with 



(A9) 



and 



rj]|i NS — Resip 



(A10) 



(All) 



As discussed in Section IA 9l we usually drop the linearized 
matter term 10ir'i' i f 2 p in equation (| AlOp above. This 
does not affect the final solution, and while removing this 
term slowed down the solution of the individual equation, 
it did improve the convergence of the overall iteration 
scheme. 

In the above equations the background scalar curva- 
ture R is given by 



R o Q BH g2 
It = o . 

r 2 
'BH 



(A12) 



and Aij is computed from equation 1)3.26)1 (see Section 
El below). 



4. The lapse equation 

It is convenient to introduce a new variable 

N = aip, (A13) 
in terms of which equation 13.25fl can be rewritten as 

D 2 N - -N*p- S A 2 - —N^K 2 - -NR - ^(3' l D t k 
8 12 8 

= 2TrN^ 4 f 2 e (p + 2S). 



(A14) 



Similar to the split (| A3|) we write N as contributions 
from the black hole and the neutron star 



N = N BB + N NS . 



(A15) 



In the absence of the neutron star TVbh reduces to oibh, 
so we may identify 



A^bh = «BH 



(A16) 
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The split (|Al"5)) can now be inserted into (|A14(i . which 
results in an equation of the form (|A1|) with 

u Nns = -^- 8 i 2 - ^K 2 -\R- 2^fl{p + 25) 

(A17) 

and 

+ la BH R + i^fi l b l K + 2na lm ip i f 2 e {p + 2S) 



where 



h W BB H 2 

U QBH = 72 • 

r BH 



(A19) 



The vector Laplacian of [3 BH is given by 



, 8a BB H 

a lPbh 7na — 

°'bh 



(2 + 12H + 12H 2 )P 



(A26) 



We note that in the absence of a neutron star, TVns = 
■0ns = SI = f — 0, equation (|A22|I is satisfied analyti- 
cally with (3^ s = 0. 

The next goal is to bring the operator 



1 



NS 



(A27) 



into the form (|A1|) . The mixed second derivatives can be 
eliminated by writing /3^ s as the sum of a vector and the 
gradient of a scalar 



/3ns = W* + l -D l U. 



(A28) 



5. The Momentum Constraint 

We first split the shift vector into contributions from 
the black hole, the neutron star, and a rotational term 



where 



C = e ijk Cl jXk = Q(-y,x- x rot ,0). 



(A20) 



(A21) 



In equation (|3.24(1 . the analytic contributions can be 
moved to the right hand side, leaving an equation for 
the neutron star contribution 

A L ^ S = 2A ij Dj(a%l)~ 6 ) + ^aD l K 

+ 167ra^f 2 e f - A L [3 BB - A L e 
ee S\ (A22) 

This equation can be solved with appropriate boundary 
conditions for /3^ s alone. 

The individual terms in l|A22|) are computed as follows. 
The first term on the right hand side is written as 

2A l3 bj{a^- 6 ) = 2A tj b J {Nijr 7 ) (A23) 
= 2A l3 -ip- 7 {b 3 a BU + DjN m ~ 7Ni(;~' 1 bjip). 

Here the first term can be computed analytically from 
the Kerr-Schild background 



f) n - Q BH# , 

UjCtBH — — tj 

r B H 



The gradient of K in (|A22|I is written as 



4 ... 47V - - . 

aD*K = —D l K 
3 3tp 



a B H + N m 8a BB H 



(A24) 



(A25) 



(2 + 10H + 9H 2 )P 



ip 3f 



With this choice, equation l|A22(l can be written as 

D 2 U = -D.W 1 (A29) 

and 



D 2 W l = S i - R\ I W j + -t)W ) . 



(A30) 



Since U is a scalar, equation l|A29|) is already in the form 
(|A1|) . and can be solved immediately with 



uu = 

r v = -DtW* 



(A31) 
(A32) 



Equation l|A30|) is more complicated. Writing out the 
covariant Laplace operator acting on a vector, we have 

D 2 W l = f j W\ ij -W\ k f k +2f j W r ^f l mj (A33) 

+w m (rt i mjti - f l mk t k + rf k mj f l kl ). 

It is convenient to define 

(A34) 



m — / mj,i mk ' > mj ki 



2a BB H 



((1 + 2H) 2 ri l m ~ (3 + 13-ff + \2H 2 )lH m ) 



' BH 

so that equation l|A33|) becomes 

D 2 W l = f j W l tj - W l k f k + 2f j W n lf l mj + W m A l m . 

(A35) 

To obtain the equation for W l , the third term and the 
pieces of the fourth term not containing W l can now be 
moved to the right hand side of equation (|A30() . The 
term K 1 i W l (no summation) on the right hand side of 
(|A30|I can also be absorbed in a u W i on the left hand 
side, which leaves us with an equation for W l in the form 
((Al"1) with 



u w , = R\ + A 1 



no summation 
1 



r W i 



R\W J - -R l ,D l U - 2W ,r ]T 



BH 



-W j A l 



U * i). 



(A36) 



(A37) 
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Hero wc have defined 



m — / rnj 
^BH 

We also note that 



(A38) 



;2(l + 2fl)lV m -(3 + 4ff)l'l i g 



*i = 



BH 



(A39) 



((1 + 6H + m 2 )rf j - (3 + lttff + 8H 2 )l l lj 



In S\ A L £ { can be computed from (|3.8|l and HA33|> . 
The partial derivatives can be computed trivially from 
(|A21|) ; in particular we have ^ = 0. Because of the 

symmetries we also have = ®i which simplifies the 

calculation. For extreme mass ratios, when a; rot = XbHj 
Ai£ l = identically. Equation l|A30|l can now be solved 
iteratively for the components W % . 



6. The Extrinsic Curvature 

With a new solution for the shift [3 Z , the extrinsic cur- 
vature A i: > can be computed from equation (|3.26|) . In 
this expression, the shift decomposition l|A20(> has to be 
inserted. The background contribution of /3g H can be 
computed analytically 



(Lp, 



BH 



(A40) 

4Q BH g (( 2 + 7H + 6H 2 )^ - (6 + 13H + 6H 2 )PP) . 
3r B H 

The contributions of the rotational term £ l are most eas- 
ily added by calculating the partial derivatives directly 
from l|A2ip . and computing the connection terms to- 
gether with those for the neutron star terms /3^ g . 



7. The integrated Euler equation 

As in the Newtonian case we first evaluate the inte- 
grated Euler equation (|3.2U|) at the three points xa, xb 
and xc to find the three eigenvalues C, f2 and f e . All 
quantities in the integrated Euler equation have to be ex- 
pressed in terms of our decompositions l|A3(l . (|A15(1 and 
(|A20|) . Of the three constants, only C appears directly in 
(|3.2U|> . The angular velocity Ct has also been introduced 
explicitly through the rotational shift term in l|A20(l . 
The scaling parameter r e again only appears through the 
implicit scaling of the gravitational potentials. In the 
Newtonian case, the Poisson equation (|2.15|) immediately 
implied a scaling for 0ns- The corresponding relativis- 
tic field equations, namely the Hamiltonian constraint 
(|3.23l) and the lapse equation (|3.25() , do not provide such 
a simple scaling because the non-linear terms. However, 
we can use the Newtonian limit to guess an appropriate 
scaling (compare 0.) 



In the Newtonian limit, the lapse becomes 

a-e^-l + ^l + ^Ns+^BH- (A41) 

Since the background describes a black hole, the confor- 
mal factor only includes a neutron star contributions, so 
that in the Newtonian limit 



>NS- 



(A42) 



Comparing with (|A3|) we note that in the Newtonian 
limit we have 



and from (|A1 5|> 



iV NS = aip - a B H ~ 2^ NS 



(A43) 



(A44) 



In analogy with the Newtonian problem, this suggests 
the definition 



tpNS = ipTsss/r 2 e 



NS 



Nm/f 2 e 



(A45) 



In the integrated Euler equation, ?/>ns & n d N^s can then 
be replaced with f^NS an d ^e-^NSj which explicitcly in- 
troduces the scaling parameter f e into the equation. 

Once the three eigenvalues C, Q and f e have been 
found, -0ns and -/Vns have to be rescaled to reflect the 
new value of f e . The integrated Euler equation can then 
be solved everywhere for the new density distribution q. 



8. Boundary Conditions 

In general, the construction of BHNS binaries requires 
three different kinds of boundary conditions: exterior 
boundary conditions at large separation from the binary, 
inner boundary conditions on the black hole's excision 
surface, and, having adopted equatorial symmetry (see 
Section IA 10 . a symmetry boundary condition on the 
z = plane. The symmetry conditions on the z — 
are straight-forward to implement (see Table ItTT|) . As ex- 
plained in Section Hvl the assumption of extreme mass 
ratios Mbh 3> Mns allows us to restrict the numerical 
grid to a region around the neutron star. This elimi- 
nates the need for interior excision boundary conditions, 
but it also requires imposing boundary conditions in the 
potentially strong-field regime close to the black hole. 

Given that we expect the neutron star to affect the 
spacetime only in a small neighborhood, it is reasonable 
to assume 1 fr Robin fall-off conditions for the conformal 
factor t/>ns an( i lapse -/Vns. 

More subtle arc the boundary conditions on the shift 
quantities W % and U. In a similar decomposition was 
used for the construction of binary neutron stars, and 
boundary conditions for these fields were identified from 
the asymptotic behavior of a multipole solution to a gen- 
eral, but flat vector Laplacian 39] . The construction 
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z — u 




V'NS 


SYM 


1 






N N S 


SYM 


1 






T 


U 


SYM 








T 


w x 


SYM 


1 






r 


w y 


SYM 


1 






r 


w z 


ANTI 


z 







TABLE III: Boundary conditions for the outer boundaries 
(OB) and on the z = symmetry plane (where SYM denotes 
symmetric and ANTI antisymmetric boundary conditions). 
Note the r is the separation from the origin of the coordinate 
system at the center of the neutron star. 



A ou t 


Grid size 


ibh/Mbh 


r e 


ymax 




2 


(32 x 32 x 16) 


17.82 


1.113 


0.02349 


0.01330 


3 


(48 x 48 x 24) 


17.81 


1.113 


0.02349 


0.01332 


4 


(64 x 64 x 32) 


17.81 


1.113 


0.02349 


0.01333 


6 


(96 x 96 x 48) 


17.81 


1.113 


0.02349 


0.01336 



TABLE IV: Numerical results for different locations of the 
outer boundary at constant grid resolution. We show values 
for a Mbh = 0.5, Mns = 0.05 binary at a separation of 4bh = 
—8 (see also the top line in Table 0. 

also assumed that the net momentum contained in the 
numerical grid vanishes (so that the non-vanishing an- 
gular momentum gives rise to the asymptotic fall-off.) 
Here, however, we solve a vector Laplacian that is not 
flat, we adopt boundary conditions which are not in an 
asymptotic region, and the grid containing only the neu- 
tron star contains a non-vanishing linear momentum. It 
is not clear a priori how appropriate boundary conditions 
on W l and U should be constructed in this case. We have 
experimented with several different conditions, and did 
notice some effect on the convergence properties of the 



iteration scheme. We found reasonable results with the 
conditions tabulated in Table 11111 which are motivated 
by the construction in [3j| but do not assume a vanish- 
ing linear momentum. In Table HVI we tabulate numerical 
values obtained for different locations of the outer bound- 
aries, which demonstrates that our results depend only 
very weakly on where the outer boundaries are imposed 
(as long as they are imposed sufficiently far away away 
from the neutron star and not too close to the black hole 
singularity). Once the assumption of extreme mass ra- 
tios is relaxed and the computational grid encompasses 
both the neutron star and the black hole, the boundary 
conditions of 0] can again be applied. 

9. Iteration Scheme 

The iteration scheme used for the construction of rela- 
tivistic BHNS binaries is very similar to that for Newto- 
nian binaries described in Section III El However, instead 
of having to solve only one Poisson equation, we now have 
to solve several Poisson-like equations together with the 
integrated Euler equation. These different equations can 
be solved in various different sequences. We found best 
convergence when solving the integrated Euler equation 
after each elliptic equation. 

We also found that the convergence improved by using 
underrelaxation, i.e. by using a linear combination 

T +1 = Xf n+1 + (1 - X)f\ (A46) 

where /" is any one of the functions N-^s, ?Ms > W l and 
U after n iteration steps, and f n+1 the function after 
having solved the corresponding elliptic equation. The 
parameter x is the relaxation parameter, and values of 
X < 1 define underrelaxation. We have found best results 
for x between 0.6 and 0.8. 

While the Hamiltonian constraint alone converged 
faster when solved as described in Section IA 31 we also 
found improved convergence of the overall iteration when 
we did not include the linearized matter term in equation 
(|A10|I . The effect of not including this term is similar to 
that of using underrelaxation, and reduced the coupling 
between the matter and conformal factor, which other- 
wise led to unstable behavior in some situations. 
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